clear; close all; clc;

%% 1. 模型参数
alpha = 6;       % (α, σ, μ) = (6, -0.1, 0.001)
sigma = -0.1;
mu    = 0.001;

% 耦合方程中 phi 的更新系数
epsilon = 1;     % ε = 1

% 迭代设置
N_total = 50000; % 总迭代次数
N_cut   = 10000; % 去除瞬态的步数

%% 2. (a) 固定 φ₀ = -0.25，扫描 k
phi0_fixed = -0.25;
k_range = linspace(-0.4, 0.4, 101);  % 扫描范围
E_a = zeros(1, length(k_range));

for i = 1:length(k_range)
    k_val = k_range(i);
    % 调用 Rulkov 函数，返回 [X1_save, X2_save, E]
    [~, ~, E_a(i)] = Rulkov(alpha, sigma, mu, k_val, phi0_fixed, N_total, epsilon, N_cut);
end

%% 3. (b) 固定 k = 0.1，扫描 φ₀
k_fixed = 0.1;
phi0_range = linspace(-2, 2, 101);  % 扫描范围
E_b = zeros(1, length(phi0_range));

for i = 1:length(phi0_range)
    phi_val = phi0_range(i);
    [~, ~, E_b(i)] = Rulkov(alpha, sigma, mu, k_fixed, phi_val, N_total, epsilon, N_cut);
end

%% 4. 绘图
figure;

% (a) 子图：E 随 k 变化（固定 φ₀ = -0.25）
subplot(1,2,1);
plot(k_range, E_a, 'LineWidth',1.5, 'Color',[0.7,0.3,0.0]);
xlabel('k');
ylabel('E');
title('(a) \phi_0 = -0.25');

% (b) 子图：E 随 φ₀ 变化（固定 k = 0.1）
subplot(1,2,2);
plot(phi0_range, E_b, 'LineWidth',1.5, 'Color',[0.0,0.6,0.0]);
xlabel('\phi_0');
ylabel('E');
title('(b) k = 0.1');
